1 Setup

1.1 Load packages

library("knitr")       # for knitting stuff
library("kableExtra")  # for markdown tables
library("lme4")        # for linear mixed effects models
library("broom.mixed") # for tidying mixed models results 
library("brms")        # Bayesian regression with Stan
library("corrr")       # for tidy correlation matrix
library("xtable")      # for latex tables
library("jsonlite")    # for reading json files
library("janitor")     # cleans stuff
library("Hmisc")       # bootstrapped confidence intervals
library("tidybayes")   # for Bayesian data analysis
library("png")         # adding pngs to images
library("grid")        # functions for dealing with images 
library("plotly")      # 3D scatter plot 
library("egg")         # for geom_custom
library("clValid")     # for validating clustering solutions 
# library("ggtern")      # for ternary plots
library("patchwork")   # for figure panels
library("tidyverse")   # for wrangling, plotting, etc. 

1.2 Helper functions

# setting some knitr options
opts_chunk$set(comment = "",
               fig.show = "hold")

# set the ggplot theme
theme_set(theme_classic())

# suppress summarise() grouping warning 
options(dplyr.summarise.inform = F)

# function for printing out html or latex tables 
print_table = function(data, format = "html", digits = 2){
  if(format == "html"){
    data %>% 
      kable(digits = digits) %>% 
      kable_styling()
  }else if(format == "latex"){
    data %>% 
      xtable(digits = digits) %>%
      print(include.rownames = F,
            booktabs = T)
  }
}

# root mean squared error
rmse = function(x, y){
  return(sqrt(mean((x - y)^2)))
}

actual_counterfactual_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = clipindex,
                            y = value,
                            fill = colorindex,
                            group = model)) +
    stat_summary(geom = "bar",
                 fun = "mean",
                 color = "black",
                 position = position_dodge(0.7),
                 width = 0.7) +
    stat_summary(geom = "linerange",
                 fun.data = "mean_cl_boot",
                 position = position_dodge(0.7),
                 size = 1) +
    geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.2) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    facet_grid(index_actual ~ index_counterfactual) +
    geom_text(data = df.text %>%
                drop_na(label),
              mapping = aes(x = x, y = y, label = as.character(label)),
              size = 6,
              position = position_dodge(width = .7),
              hjust = 0.5) +
    coord_cartesian(xlim = c(0.5, 2.5),
                    ylim = c(-5, 105),
                    expand = F,
                    clip = "off") +
    labs(y = ylabel) +
    theme(text = element_text(size = 20),
          legend.position = "none",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
          panel.spacing.y = unit(0.8, "cm"),
          plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
          axis.title.x = element_blank(),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  # print(p)
}

actual_counterfactual_threeballs_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = ball,
                            y = value,
                            group = model,
                            fill = colorindex)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = "black",
                 position = position_dodge(0.9),
                 width = 0.9) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 color = "black",
                 position = position_dodge(0.9)) +
    geom_point(data = x %>% 
                 mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(jitter.width = 0.3,
                                               jitter.height = 0,
                                               dodge.width = 0.9),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.15) +
    facet_wrap(~clip, ncol = 8) +
    geom_text(data = df.text, 
              mapping = aes(x = x,
                            y = y,
                            label = as.character(label)),
              size = 5, 
              position = position_dodge(width = .9), 
              hjust = 0.5, 
              fontface = "bold") +
    coord_cartesian(xlim = c(0.5, 2.5), 
                    ylim = c(-5, 120), 
                    expand = F,
                    clip = "off") +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    scale_color_manual(values = c("red2", "green2", "black")) +
    labs(y = ylabel) +
    theme(text = element_text(size = 16),
          legend.position = "none",
          axis.title.x = element_blank(),
          panel.spacing.y = unit(0.3, "cm"),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
}

2 Experiment 1: 2 balls

2.1 Clips

2.2 Read in data

# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv", 
                          header = T,
                          stringsAsFactors = F)

# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# Information about each clip 
df.exp1.clipinfo = tibble(clip = 1:18,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 
                                             1, 1, 0, 1, 1, 1, 1, 1, 1),
                          outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1, 
                                                     0, 1, 1, 0, 0, 1, 0, 1, 1),
                          index_actual = rep(c("actual miss", 
                                               "actual close", 
                                               "actual hit"), each = 6),
                          index_counterfactual = rep(rep(c("counterfactual miss",
                                                           "counterfactual close", 
                                                           "counterfactual hit"),
                                                         each = 2),
                                                     3)) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp1.causal.long = df.exp1.causal$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

2.2.1 Demographics

# counterfactual condition 
df.exp1.counterfactual %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
41 32 8.8 19 7.5 3.12
# causal condition 
df.exp1.causal %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
41 33 10.3 21 9.65 6.58

2.3 Read in model predictions

2.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]

df.exp1.model = files %>%
  map_dfr(~ read.csv(str_c(path, .))) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp1.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
  group_by(clip) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>%
  ungroup() %>%
  left_join(df.exp1.clipinfo, by = "clip")

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp1.counterfactual.model = df.exp1.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, 
                   ~ sum((.x$prediction*100 - 
                            df.exp1.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

2.3.2 Causal condition

df.exp1.causal.model = df.exp1.counterfactual.long %>% 
  group_by(clip) %>% 
  summarize(empirical = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.counterfactual.model %>% 
              select(-c(noise, sse)),
            by = "clip") %>% 
  rename(simulation = prediction) %>% 
  mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
         simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))

2.4 Counterfactual condition

2.4.1 Plots

set.seed(1)
df.plot = df.exp1.counterfactual.long %>% 
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.)/2)) %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", "actual close", "actual hit")),
         index_counterfactual = factor(index_counterfactual,
                                       levels = c("counterfactual miss",
                                                  "counterfactual close",
                                                  "counterfactual hit")),
         model = factor(model, levels = c("rating", "prediction"))) %>% 
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12


df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(x = df.plot, ylabel = "counterfactual judgment")
print(p)

# ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

2.4.2 Stats

2.4.2.1 Model evaluation

df.exp1.counterfactual.means %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise), 
            simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  print_table()
noise simulation_r simulation_rmse
1 0.97 10.95

2.5 Causal condition

2.5.1 Plots

set.seed(1)

model.name = c("empirical")
# model.name = c("simulation")

df.plot = df.exp1.causal.long %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp1.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, simulation, empirical),
               names_to = "model",
               values_to = "value") %>% 
  filter(model %in% c(model.name, "rating")) %>%
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(
    colorindex = ifelse(outcome_actual == 1, 2, 1),
    colorindex = ifelse(model != "rating", 3, colorindex),
    colorindex = as.factor(colorindex),
    index.actual = factor(index_actual, 
                          levels = c("actual miss", "actual close", "actual hit")),
    index.counterfactual = factor(index_counterfactual, 
                                  levels = c("counterfactual miss", 
                                             "counterfactual close", 
                                             "counterfactual hit"))) %>%
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) # swap clips 11 and 12

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(df.plot, "causal judgment")
print(p)

# ggsave("../../figures/plots/exp1_causal_bars.pdf",
#        width = 8,
#        height = 6)

2.5.2 Stats

2.5.2.1 Model evaluation

df.exp1.causal.long %>% 
  mutate(rating = abs(rating)) %>% 
  group_by(clip,
           outcome_actual, 
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>% 
  summarize(empirical_r = cor(rating, empirical),
            empirical_rmse = rmse(rating, empirical),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>% 
  print_table()
empirical_r empirical_rmse simulation_r simulation_rmse
0.96 8.57 0.93 15.15

2.5.2.2 Bayesian regression

# random intercepts model 
fit.brm_exp1_causal = brm(formula = rating ~ 1 + index_actual + index_counterfactual + 
                            outcome_actual + (1 + index_actual + index_counterfactual + 
                                                outcome_actual | participant),
                          data = df.exp1.causal.long,
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp1_causal")

2.5.2.3 Specific hypotheses

df.exp1.causal.posterior = df.exp1.causal.long %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp1_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)

func_posterior_difference = function(data, clip1, clip2){
  data %>% 
    mutate(difference = .data[[clip1]] - .data[[clip2]]) %>% 
    pull(difference) %>% 
    mean_hdci() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>% 
    print()  
}

func_posterior_difference(data = df.exp1.causal.posterior,
                          clip1 = "8",
                          clip2 = "13")
            difference
1 (0.26 [-6.51, 7.31])

2.5.3 Tables

fit.brm_exp1_causal %>% 
  tidy(effects = "fixed") %>% 
  mutate(term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = conf.low,
         `upper 95% HDI` = conf.high) %>% 
  # print_table(format = "latex")
  print_table()
Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
effect component term estimate std.error lower 95% HDI upper 95% HDI
fixed cond (intercept) -18.33 3.36 -24.86 -11.66
fixed cond index_actualactualclose 4.31 3.40 -2.54 10.94
fixed cond index_actualactualhit 4.06 4.84 -5.51 13.27
fixed cond index_counterfactualcounterfactualclose -17.52 3.37 -24.36 -10.84
fixed cond index_counterfactualcounterfactualhit -61.37 4.38 -69.95 -52.89
fixed cond outcome_actual 92.36 5.11 82.21 102.38

3 Experiment 2: Brick & Teleport

3.1 Clips

3.2 Read in data

df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
                                header = T,
                                stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
                                        header = T,
                                        stringsAsFactors = F)

# Information about each clip

df.exp2.clipinfo = tibble(clip = 1:20,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 
                                             0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
                          outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 
                                                     1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
                          index_actual = c(rep(c("actual miss", 
                                                 "actual close",
                                                 "actual hit"),
                                               each = 6),
                                           rep("practice", 2)),
                          index_counterfactual = c(rep(rep(c("counterfactual miss",
                                                             "counterfactual close",
                                                             "counterfactual hit"),
                                                           each = 2), 
                                                       3),
                                                   rep("practice", 2))) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp2.causal_first.long = df.exp2.causal_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("causal", "counterfactual"), each = 20),
                        max(participant)),
         condition = "causal_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("counterfactual", "causal"), each = 20),
                        max(participant)),
         condition = "counterfactual_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
  bind_rows(df.exp2.counterfactual_first.long %>% 
              mutate(participant = participant + 
                       max(df.exp2.causal_first.long$participant)))

3.2.1 Demographics

# counterfactual condition 
df.exp2.counterfactual_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time) %>% 
  bind_rows(df.exp2.causal_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time)) %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
82 35 12.2 45 21.25 5.11

3.3 Read in model predictions

3.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]

df.exp2.model = files %>%
  map_dfr(~ read.csv(str_c(path, .))) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp2.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp2.counterfactual.means = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  # mutate(rating = abs(rating)) %>% 
  group_by(clip) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.clipinfo, by = "clip")

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp2.counterfactual.model = df.exp2.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, 
                   ~ sum((.x$prediction*100 - df.exp2.counterfactual.means$rating) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

3.3.2 Causal condition

# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
  filter(question == "counterfactual",
    clip <= 18) %>%
  group_by(clip, condition) %>%
  summarize(rating = mean(rating)) %>%
  pivot_wider(names_from = condition,
              values_from = rating) %>% 
  left_join(df.exp2.combined.long %>%
              filter(question == "counterfactual", clip <= 18) %>%
              group_by(clip) %>%
              summarize(combined = mean(rating)),
            by = "clip") %>% 
  left_join(df.exp2.counterfactual.model,
            by = "clip") %>% 
  select(-c(sse, noise)) %>% 
  rename(simulation = prediction) %>% 
  mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined),
            .funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))

3.4 Counterfactual condition

3.4.1 Plots

set.seed(1)

df.plot = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.counterfactual.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")),
         model = factor(model,
                             levels = c("rating", "prediction")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"),
                     each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(df.plot, "counterfactual judgment")
print(p)

# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

3.4.2 Stats

3.4.2.1 Model evaluation

df.exp2.counterfactual.means %>% 
  left_join(df.exp2.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise),
            simulation_r = cor(rating, prediction),
            simulation_rmse = rmse(rating, prediction)) %>% 
  print_table()
noise simulation_r simulation_rmse
0.9 0.96 13.46

3.4.2.2 Bayesian regressions

3.4.2.2.1 Check for potential order effects
# random intercepts & slopes model 
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition * 
                                    (index_actual + index_counterfactual) + 
                                    (1 + index_actual + 
                                       index_counterfactual | participant),
                                  data = df.exp2.combined.long %>% 
                                    na.omit() %>% 
                                    filter(question == "counterfactual"),
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp2_counterfactual")

3.4.3 Tables

fit.brm_exp2_counterfactual %>% 
  tidy(effects = "fixed") %>% 
  mutate(term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = conf.low,
         `upper 95% HDI` = conf.high) %>% 
  # print_table(format = "latex")
  print_table()
Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
effect component term estimate std.error lower 95% HDI upper 95% HDI
fixed cond (intercept) 13.38 2.58 8.30 18.38
fixed cond conditioncounterfactual_first -1.11 3.65 -7.99 6.04
fixed cond index_actualactualclose -1.48 2.78 -6.79 3.95
fixed cond index_actualactualhit -4.42 2.40 -9.15 0.20
fixed cond index_counterfactualcounterfactualclose 37.94 3.33 31.43 44.60
fixed cond index_counterfactualcounterfactualhit 73.80 4.17 65.54 81.99
fixed cond conditioncounterfactual_first:index_actualactualclose -0.55 3.84 -8.29 6.79
fixed cond conditioncounterfactual_first:index_actualactualhit -0.85 3.41 -7.60 5.87
fixed cond conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.23 4.79 -8.37 10.75
fixed cond conditioncounterfactual_first:index_counterfactualcounterfactualhit -0.02 5.92 -11.96 11.59

3.5 Causal condition

3.5.1 Plots

set.seed(1)

func_exp2_causal_plot = function(condition.name, model.name){

df.plot = df.exp2.combined.long %>%
  filter(question == "causal", clip <= 18) %>%
  filter(condition %in% condition.name) %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c("rating", all_of(model.name)),
               names_to = "model",
               values_to = "value") %>% 
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual, 
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"),
                     each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(df.plot, "causal judgment")
}

plot.list = list(func_exp2_causal_plot(condition.name = "counterfactual_first",
                                       model.name = "combined"),
                 func_exp2_causal_plot(condition.name = "causal_first",
                                       model.name = "combined"))

wrap_plots(plot.list, ncol = 2)

# condition.name = "counterfactual_first"
# condition.name = "causal_first"
# model.name = "combined"
# 
# func_exp2_causal_plot(condition.name = condition.name,
#                       model.name = model.name)

# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
#        width = 8,
#        height = 6)
Left panel: Counterfactual judgments first; right panel: Causal judgments first

Figure 3.1: Left panel: Counterfactual judgments first; right panel: Causal judgments first

3.5.2 Stats

3.5.2.1 Model evaluation

df.data = df.exp2.combined.long %>%
  filter(question == "causal",
         clip <= 18) %>%
  mutate(rating = abs(rating)) %>%
  group_by(condition,
           clip,
           outcome_actual,
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual"))

df.data %>% 
  group_by(condition) %>% 
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined)) %>% 
  print_table()
condition empirical_r empirical_rmse
causal_first 0.92 13.78
counterfactual_first 0.99 5.27
df.data %>%
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>%
  print_table()
empirical_r empirical_rmse simulation_r simulation_rmse
0.95 10.43 0.92 19.09

3.5.2.2 Bayesian regression

# random intercepts + slopes model 
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition * 
                            (index_actual + index_counterfactual + outcome_actual) +
                            (1 + index_actual + index_counterfactual + 
                               outcome_actual | participant),
                          data = df.exp2.combined.long %>% 
                            filter(question == "causal",
                                   clip <= 18),
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp2_causal")

3.5.2.3 Specific hypotheses

df.exp2.causal.posterior = df.exp2.combined.long %>% 
  filter(question == "causal",
         clip <= 18) %>% 
  group_by(clip, index_actual, index_counterfactual,
           outcome_actual, condition) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp2_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, condition, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)

# difference between clips 13 and 14 versus 8 
df.exp2.causal.posterior %>% 
    mutate(difference = (`13`  + `14`) / 2 - `8`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))
               difference
1 (-0.55 [-12.09, 10.97])
# difference between clips 5 and 6 versus 11 
df.exp2.causal.posterior %>% 
    mutate(difference = (`5`  + `6`) / 2 - `11`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))
             difference
1 (-5.48 [-13.8, 2.73])

3.5.3 Tables

3.5.3.1 Bayesian regression

fit.brm_exp2_causal %>% 
  tidy(effects = "fixed") %>% 
  mutate(term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = conf.low,
         `upper 95% HDI` = conf.high) %>% 
  # print_table(format = "latex")
  print_table()
Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
effect component term estimate std.error lower 95% HDI upper 95% HDI
fixed cond (intercept) -24.07 3.50 -30.87 -17.24
fixed cond conditioncounterfactual_first 11.47 4.91 1.62 21.22
fixed cond index_actualactualclose 8.17 3.44 1.47 14.79
fixed cond index_actualactualhit 13.42 4.90 3.74 22.69
fixed cond index_counterfactualcounterfactualclose -30.20 3.72 -37.61 -23.19
fixed cond index_counterfactualcounterfactualhit -50.29 4.60 -59.57 -41.47
fixed cond outcome_actual 98.53 5.40 87.81 108.98
fixed cond conditioncounterfactual_first:index_actualactualclose -5.39 4.94 -15.00 4.35
fixed cond conditioncounterfactual_first:index_actualactualhit -16.97 7.02 -30.40 -3.33
fixed cond conditioncounterfactual_first:index_counterfactualcounterfactualclose -8.13 5.17 -18.18 2.17
fixed cond conditioncounterfactual_first:index_counterfactualcounterfactualhit -25.51 6.44 -38.16 -12.98
fixed cond conditioncounterfactual_first:outcome_actual 10.72 7.60 -4.16 26.02

4 Experiment 3: 3 balls

4.1 Clips

4.2 Read in data

# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
                          outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 
                                           0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
                          outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
                          outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
                          outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                                           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                          clipindex = rep(1:2, 16))

# COUNTERFACTUAL JUDGMENTS 
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty",
                      "smoothness", "time", "condition"),
                    n()/6),
         participant = rep(1:(n()/6), each = 6)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, condition, gender, age, time, smoothness, difficulty)

# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "rating"), n()/2),
         participant = rep(1:(n()/64), each = 64),
         order = rep(rep(1:32, each = 2), n()/64)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  left_join(df.exp3.counterfactual.demographics %>% 
              select(participant, ball = condition),
            by = "participant") %>% 
  select(participant, ball, clip, everything()) %>% 
  arrange(participant, clip)

# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
                          header = T,
                          stringsAsFactors = F)

# demographics 
df.exp3.causal.demographics = df.exp3.causal$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", 
                      "counterbalance", "replayType", "experiment"),
                    n()/8),
         participant = rep(1:(n()/8), each = 8)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, gender, age, counterbalance, time, smoothness, difficulty)

# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
         participant = rep(1:(n()/(32*6)), each = (32*6)),
         order = rep(rep(1:32, each = 6), n()/(32*6))) %>% 
  filter(!name %in% c("x", "y", "z")) %>% 
pivot_wider(names_from = name,
            values_from = value) %>% 
  left_join(df.exp3.causal.demographics %>% 
              select(participant, counterbalance),
            by = "participant") %>% 
  mutate(A = ifelse(counterbalance == 1, rating1, rating2),
         B = ifelse(counterbalance == 1, rating2, rating1)) %>% 
  select(-rating1, -rating2) %>% 
  pivot_longer(cols = c(A, B),
               names_to = "ball",
               values_to = "rating") %>% 
  select(participant, clip, ball, rating, order) %>% 
  arrange(participant, clip, ball)

4.2.1 Demographics

# counterfactual condition 
df.exp3.counterfactual.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
80 33 10.1 34 18.08 4.63
# causal condition 
df.exp3.causal.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
41 34 10.5 21 21.19 4.96

4.3 Read in model predictions

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]

df.exp3.model = files %>%
  map_dfr(~ read_csv(str_c(path, .))) %>% 
  rename(clip = trial)

4.3.1 Counterfactual condition

# calculate mean counterfactual judgments 
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

# find noisy simulation model that best predicts the mean counterfactual 
# judgments by minimizing the sum of squared errors 
df.exp3.model.counterfactual = df.exp3.model %>% 
  group_by(noise) %>% 
  select(clip, contains("whether"), noise) %>% 
  pivot_longer(cols = c(A_whether, B_whether),
               names_to = "ball",
               values_to = "prediction") %>% 
  mutate(ball = str_remove(ball, "_whether")) %>% 
  arrange(noise, clip, ball) %>% 
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data,
                   ~ sum((.x$prediction*100 -
                            df.exp3.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

4.3.2 Causal condition

# calculate mean causal judgments 
df.exp3.causal.means = df.exp3.causal.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

df.exp3.model.causal = df.exp3.model %>% 
  # take into account difference-making 
  mutate_at(.vars = vars(A_how:A_robust),
            .funs = list(~ . * A_difference)) %>% 
  mutate_at(.vars = vars(B_how:B_robust),
            .funs = list(~ . * B_difference)) %>% 
  # choose model based on best fit with counterfactual data 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>% 
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(clip, ball:robust)

4.3.3 Heuristic model

l.features = fromJSON("data/features.json")

df.features.initial = l.features[["appearance"]] %>%
  cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
  setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
             "ballA_velx",
             "ballA_vely",
             "ballB_velx",
             "ballB_vely",
             "ballE_velx",
             "ballE_vely")) %>%
  mutate(clip = 1:nrow(.)) %>%
  pivot_longer(cols = starts_with("ball")) %>%
  separate(name, into = c("ball", "property")) %>%
  pivot_wider(names_from = property, values_from = value) %>% 
  mutate(ball = str_remove_all(ball, pattern = "ball"))

# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")

for (i in 1:32) {
  df.tmp = data.frame()
  for (j in 1:length(ballnames)) {
    ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
    if (ncollisions > 0) {
      tmp = data.frame(
        ball = ballnames[j],
        object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
        time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
        pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
        pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
        post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
        post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
        ncollision = 1:ncollisions
      )
      df.tmp = df.tmp %>%
        rbind(tmp)
    }
  }
  df.features.collisions = df.features.collisions %>%
    rbind(df.tmp %>%
            mutate(clip = i) %>%
            select(clip, ball, everything()))
}

# find end of clip
tmp = df.features.collisions %>%
  group_by(clip) %>%
  filter(ball == "ballE",
         str_detect(object, "goal|Left|Right")) %>%
  group_by(clip) %>%
  mutate(endclip = max(time)) %>%
  select(clip, endclip) %>%
  ungroup() %>%
  distinct()

# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
  left_join(tmp, by = "clip") %>%
  mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
  group_by(clip) %>%
  filter(time <= endclip)

# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
  filter(ball == "E") %>%
  mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
  select(clip, E.moving)

# contact with ball E
tmp.contactE = df.features.collisions %>%
  filter(object == "ballE") %>%
  mutate(ball = factor(ball,
                       levels = c("ballA", "ballB"),
                       labels = c("A", "B"))) %>%
  mutate(contactE = 1) %>%
  select(clip, ball, contactE) %>%
  group_by(clip, ball) %>%
  summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate(contactE = ifelse(is.na(contactE), 0, contactE))

# number of collisions
tmp.ncollisions = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  mutate(ball = factor(ball,
                       levels = c("ballA", "ballB"),
                       labels = c("A", "B"))) %>%
  group_by(clip, ball) %>%
  count() %>%
  rename(ncollision = n)

# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  count(clip, ball, object) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
  group_by(clip) %>%
  summarize(AE = any(AE %>% as.logical()) * 1,
            BE = any(BE %>% as.logical()) * 1) %>%
  mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
         B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
  select(clip, A, B) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "exclusive")

# impact
func_angle = function(x, y) {
  dot.prod = x %*% y
  norm.x = norm(x, type = "2")
  norm.y = norm(y, type = "2")
  theta = acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}

# impact on ball E 
tmp.impactE = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ball"),
         ball == "ballE") %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely),
                                     c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely,
                                           post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball,
                       levels = c("ballA", "ballB"),
                       labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(E.speed.diff = speed.diff,
         E.direction.diff = direction.diff)

# total impact 
tmp.impactTotal = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ballA|ballB")) %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely),
                                     c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely,
                                           post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball,
                       levels = c("ballA", "ballB"),
                       labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(total.speed.diff = speed.diff,
         total.direction.diff = direction.diff)

# transfer of force 
tmp.transfer = df.features.collisions %>%
  filter(str_detect(ball, "ballA|ballB"),
         str_detect(object, "ball")) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
         AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
  mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
  filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
  arrange(clip, time) %>%
  group_by(clip) %>%
  summarize(A = any(!is.na(AE)),
            B = any(!is.na(BE)),
            A = ifelse(any(!is.na(AB)) &
                         any(!is.na(BE)) &
                         max(AB, na.rm = T) < min(BE, na.rm = T),
                       T,
                       A),
            B = ifelse(any(!is.na(AB)) &
                         any(!is.na(AE)) &
                         max(AB, na.rm = T) < min(AE, na.rm = T),
                       T,
                       B)) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "transfer") %>% 
  arrange(clip, ball)

# collect features
df.features = df.features.initial %>%
  filter(ball != "E") %>%
  mutate(ball = factor(ball)) %>%
  mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
         speed = abs(velx) + abs(vely)) %>%
  left_join(tmp.contactE) %>%
  left_join(tmp.impactE) %>%
  left_join(tmp.impactTotal) %>%
  left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
  left_join(tmp.movingE) %>%
  left_join(tmp.exclusive) %>%
  select(-c(appearance, velx, vely))

# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])

4.4 Counterfactual condition

4.4.1 Plots

4.4.1.1 Bar graph

set.seed(1)

df.plot = df.exp3.counterfactual.long %>%
  left_join(df.exp3.model.counterfactual %>%
              select(clip, ball, prediction) %>%
              mutate(ball = as.factor(ball)),
            by = c("clip", "ball")) %>%
  left_join(df.exp3.clipinfo,
            by = "clip") %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "prediction", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "prediction"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>%
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

p = actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
#        width = 12,
#        height = 6)

4.4.2 Stats

4.4.2.1 Model evaluation

# best fitting model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.model.counterfactual,
            by = c("clip", "ball")) %>% 
  summarize(simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  print_table()
simulation_r simulation_rmse
0.86 19.84
# deterministic model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_a, outcome_b) %>% 
              pivot_longer(cols = -clip,
                           names_to = "ball",
                           values_to = "outcome") %>% 
              mutate(ball = factor(ball,
                                   levels = c("outcome_a", "outcome_b"),
                                   labels = c("B", "A"))),
            by = c("clip", "ball")) %>% 
  mutate(outcome = outcome * 100) %>% 
  summarize(simulation_r = cor(rating_mean, outcome),
            simulation_rmse = rmse(rating_mean, outcome)) %>% 
  print_table()
simulation_r simulation_rmse
0.83 30.28

4.5 Causal condition

4.5.1 Stats

4.5.1.1 Bayesian mixed effects model

df.data = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# using whether-causation from the model 
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether + 
                              (1 + whether | participant),
                                  data = df.data,
                                  save_all_pars = T,
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp3_causal_w")

fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how +
                               (1 + whether + how | participant),
                                   data = df.data,
                                   save_all_pars = T,
                                   cores = 2,
                                   seed = 1,
                                   file = "cache/brm_exp3_causal_wh")

fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient +
                                (1 + whether + how + sufficient | participant),
                                    data = df.data,
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_whs")

4.5.1.2 Model comparison

fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
                                            criterion = c("loo", "waic"),
                                            reloo = T)
No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
                                             criterion = c("loo", "waic"),
                                             reloo = T)
No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
                                              criterion = c("loo", "waic"),
                                              reloo = T)
No problematic observations found. Returning the original 'loo' object.
loo_compare(fit.brm_exp3_causal_w,
            fit.brm_exp3_causal_wh,
            fit.brm_exp3_causal_whs)
                        elpd_diff se_diff
fit.brm_exp3_causal_whs    0.0       0.0 
fit.brm_exp3_causal_wh  -125.2      15.6 
fit.brm_exp3_causal_w   -426.2      31.1 

4.5.1.3 Heuristic model

  • z-scored predictors and invidiual responses
# Regression based on features
df.regression.features = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  left_join(df.features %>% 
              mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
            by = c("clip", "ball")) %>%
  clean_names() %>% 
  mutate(e_moving = 1 - e_moving)

# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving + 
                                      speed + 
                                      contact_e + 
                                      e_speed_diff + 
                                      e_direction_diff + 
                                      total_speed_diff + 
                                      total_direction_diff + 
                                      transfer + 
                                      e_moving + 
                                      exclusive + (1 | participant),
                                    prior = priors,
                                    data = df.regression.features %>% 
                                      select(-c(clip, ball, order, clipindex,
                                                contains("outcome"))),
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant) 
   Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     8.53      1.22     6.44    11.18 1.00     1163     1936

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept               49.79      1.46    46.91    52.74 1.00     1129
moving                   0.22      0.22     0.01     0.81 1.00     3399
speed                    2.06      0.85     0.44     3.67 1.00     2319
contact_e                0.38      0.35     0.01     1.29 1.00     3201
e_speed_diff             0.12      0.12     0.00     0.46 1.00     4435
e_direction_diff         1.04      0.74     0.04     2.71 1.00     2202
total_speed_diff         2.21      0.95     0.41     4.17 1.00     2208
total_direction_diff     3.99      0.91     2.14     5.68 1.00     2718
transfer                15.59      0.79    14.06    17.14 1.00     3555
e_moving                 0.18      0.17     0.01     0.62 1.00     3296
exclusive                4.39      0.73     2.97     5.81 1.00     3526
                     Tail_ESS
Intercept                1854
moving                   1547
speed                    1266
contact_e                2020
e_speed_diff             2268
e_direction_diff         1685
total_speed_diff         1393
total_direction_diff     1742
transfer                 2992
e_moving                 1738
exclusive                2538

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    33.21      0.46    32.33    34.10 1.00     5313     2845

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

4.5.1.4 Predictions for mean causal judgments

func_fit_data = function(fit){
  fit %>% 
    fitted(newdata = df.exp3.model.causal %>% 
             left_join(df.features %>% 
                         mutate_at(.vars = vars(moving:exclusive),
                                   .funs = ~ scale(.)[,]),
                       by = c("clip", "ball")) %>%
             clean_names() %>% 
             mutate(e_moving = 1 - e_moving),
           re_formula = NA) %>% 
    as_tibble() %>% 
    pull(Estimate)
}

df.exp3.causal.regression = df.exp3.causal.means %>% 
  mutate(w = func_fit_data(fit.brm_exp3_causal_w),
         wh = func_fit_data(fit.brm_exp3_causal_wh),
         whs = func_fit_data(fit.brm_exp3_causal_whs),
         heuristic = func_fit_data(fit.brm_exp3_causal_heuristic))

4.5.1.5 Model evaluation

prediction = "whs"

df.exp3.causal.regression %>% 
  summarize(simulation_r = cor(rating_mean, .[[prediction]]),
            simulation_rmse = rmse(rating_mean, .[[prediction]])) %>% 
  print_table()
simulation_r simulation_rmse
0.87 13.07

4.5.1.6 Individual participant regressions

df.fit = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# priors 
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline = 
  brm(formula = rating ~ 1,
      data = df.fit %>% 
        filter(participant == 1),
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_baseline"))

fit.brm_exp3_causal_individual_w = 
  brm(formula = rating ~ 1 + whether,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_w"))

fit.brm_exp3_causal_individual_wh = 
  brm(formula = rating ~ 1 + whether + how,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_wh"))

fit.brm_exp3_causal_individual_whs = 
  brm(formula = rating ~ 1 + whether + how + sufficient,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_whs"))

# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
  group_by(participant) %>%
  nest() %>%
  ungroup() %>% 
  # fit model for each participant
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.brm_exp3_causal_individual_baseline,
                                          newdata = .x)),
         fit_w = map(.x = data,
                     .f = ~ update(fit.brm_exp3_causal_individual_w,
                                   newdata = .x)),
         fit_wh = map(.x = data,
                      .f = ~ update(fit.brm_exp3_causal_individual_wh,
                                    newdata = .x)),
         fit_whs = map(.x = data,
                       .f = ~ update(fit.brm_exp3_causal_individual_whs,
                                     newdata = .x))) %>% 
  mutate(fit_baseline = map(.x = fit_baseline,
                            ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_w = map(.x = fit_w,
                     ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_wh = map(.x = fit_wh,
                      ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_whs = map(.x = fit_whs,
                       ~ add_criterion(.x, criterion = c("loo", "waic"))),
         r_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ cor(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         r_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ cor(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         r_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ cor(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         rmse_w = map2_dbl(.x = data,
                           .y = fit_w,
                           .f = ~ rmse(.x$rating, .y %>% 
                                         fitted() %>% 
                                         .[, 1])),
         rmse_wh = map2_dbl(.x = data,
                            .y = fit_wh,
                            .f = ~ rmse(.x$rating, .y %>% 
                                          fitted() %>% 
                                          .[, 1])),
         rmse_whs = map2_dbl(.x = data,
                             .y = fit_whs,
                             .f = ~ rmse(.x$rating, .y %>% 
                                           fitted() %>% 
                                           .[, 1])),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           w = fit_w,
                                           wh = fit_wh,
                                           whs = fit_whs),
                                 .f = ~ loo_compare(..1, ..2, ..3, ..4)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3", "..4"),
                             labels = c("baseline", "w", "wh", "whs")))

save(list = c("df.exp3.causal.individual_fit"),
     file = "data/exp3_causal_individual_fit.RData")
4.5.1.6.1 Load individual participant regressions
load("data/exp3_causal_individual_fit.RData")

# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>% 
  count(best_model) %>% 
  print_table()
best_model n
wh 2
whs 39
4.5.1.6.2 Model performance on individual participant level
df.exp3.causal.individual_fit %>% 
  select(participant, contains("r_"), contains("rmse_")) %>% 
  pivot_longer(cols = -participant,
               names_to = c("measure", "model"),
               names_sep = "_",
               values_to = "fit") %>% 
  group_by(model, measure) %>% 
  summarize(quantiles = quantile(fit, probs = c(0.05, 0.5, 0.95)),
            prob = c(0.05, 0.5, 0.95)) %>% 
  pivot_wider(names_from = prob,
              values_from = quantiles) %>% 
  mutate(across(.fns = ~ round(., 2))) %>% 
  print_table()
model measure 0.05 0.5 0.95
w r 0.18 0.43 0.57
w rmse 29.33 36.03 43.62
wh r 0.35 0.60 0.78
wh rmse 23.82 32.19 40.60
whs r 0.39 0.67 0.79
whs rmse 22.69 31.17 38.12
4.5.1.6.3 Clusters of individual participants’ responses
set.seed(1)

df.cluster_whs = fit.brm_exp3_causal_whs %>% 
  ranef() %>% 
  .$participant %>% 
  as_tibble() %>% 
  select(contains("Estimate")) %>% 
  set_names(tolower(str_replace_all(names(.), "Estimate.", ""))) %>% 
  mutate(participant = 1:n()) %>% 
  relocate(participant)

df.cluster_whs = df.cluster_whs %>% 
  mutate(cluster = df.cluster_whs %>%
           select(-participant) %>%
           kmeans(centers = 2) %>% 
           .$cluster)

df.cluster_whs %>% 
  print_table()
participant intercept whether how sufficient cluster
1 8.30 0.44 -5.36 -0.84 1
2 0.03 -20.89 32.67 -20.83 2
3 2.27 -8.27 26.27 -9.42 2
4 -0.66 7.98 -2.96 10.74 1
5 -0.75 2.47 -5.92 0.20 1
6 5.09 -15.71 22.42 -18.09 2
7 -0.17 1.71 -3.13 0.22 1
8 -3.55 -2.49 4.88 -0.21 2
9 -4.15 -5.31 9.38 -2.57 2
10 2.10 -9.99 3.25 -10.96 2
11 2.10 7.99 -10.81 5.96 1
12 -2.00 15.83 -25.89 18.27 1
13 3.03 4.98 0.72 2.29 1
14 3.17 4.02 -11.54 1.56 1
15 -5.07 -13.56 -3.49 -11.90 2
16 -0.89 23.98 -17.31 25.30 1
17 0.35 -4.04 2.84 -3.79 2
18 3.05 -0.15 -10.15 -2.03 1
19 -2.30 -7.53 4.17 -7.76 2
20 8.44 20.15 -11.67 17.52 1
21 -1.69 1.54 -12.89 2.37 1
22 -5.38 1.33 -3.82 2.56 1
23 -0.53 1.24 -7.19 3.29 1
24 -0.43 1.91 -0.12 2.55 1
25 2.88 12.33 -12.26 9.65 1
26 0.46 12.55 1.40 11.74 1
27 -4.46 16.00 -20.07 21.43 1
28 2.21 -8.24 -1.57 -10.35 2
29 -0.78 7.61 -4.96 7.93 1
30 -3.78 -21.34 18.53 -19.98 2
31 -6.59 7.50 -12.71 12.17 1
32 -1.70 2.50 3.37 6.52 1
33 3.29 13.81 -10.56 13.26 1
34 -4.20 -12.18 9.83 -10.61 2
35 0.41 -16.26 0.37 -18.28 2
36 -1.55 -27.69 47.62 -28.75 2
37 10.64 13.67 0.43 8.39 1
38 1.35 9.88 -7.13 9.38 1
39 -5.71 -10.20 7.58 -10.99 2
40 -6.18 -7.24 8.60 -6.13 2
41 1.97 -2.12 -1.61 -0.67 1
# cluster sizes
df.cluster_whs %>% 
  count(cluster)
# A tibble: 2 x 2
  cluster     n
    <int> <int>
1       1    25
2       2    16
4.5.1.6.3.1 Check that clustering is stable

A 2-cluster solution yields a good result according to a number of different validation measures (see here for more details on the different measures).

tmp = df.cluster_whs %>% 
  select(-participant, -cluster)

fit = clValid(obj = as.matrix(tmp),
        nClust = 2:10,
        clMethods = c("kmeans"),
        validation = c("internal", "stability"))
Warning in clValid(obj = as.matrix(tmp), nClust = 2:10, clMethods =
c("kmeans"), : rownames for data not specified, using 1:nrow(data)
fit %>% 
  summary()

Clustering Methods:
 kmeans 

Cluster sizes:
 2 3 4 5 6 7 8 9 10 

Validation Measures:
                           2       3       4       5       6       7       8       9      10
                                                                                            
kmeans APN            0.0991  0.1419  0.1438  0.1181  0.2214  0.1734  0.2094  0.1964  0.2269
       AD            20.3268 15.8118 13.3966 11.5686 11.4495 10.0034  9.6513  9.2184  8.8807
       ADM            7.3494  4.1130  4.3773  2.2086  3.7885  3.1164  3.5423  3.9327  3.9086
       FOM            8.3657  7.0101  6.0226  5.4730  5.1633  4.9791  4.8082  4.8092  4.9224
       Connectivity   4.4262 15.5893 20.1163 22.5329 25.0329 28.1413 34.4710 42.6353 44.9687
       Dunn           0.1046  0.1626  0.1687  0.2509  0.2509  0.2550  0.1815  0.2524  0.2524
       Silhouette     0.4691  0.3871  0.3868  0.3738  0.3522  0.3255  0.2988  0.2846  0.2774

Optimal Scores:

             Score  Method Clusters
APN          0.0991 kmeans 2       
AD           8.8807 kmeans 10      
ADM          2.2086 kmeans 5       
FOM          4.8082 kmeans 8       
Connectivity 4.4262 kmeans 2       
Dunn         0.2550 kmeans 7       
Silhouette   0.4691 kmeans 2       

4.5.2 Plots

4.5.2.1 Bar graph

set.seed(1)
model_index = "whs"

# model predictions 
model_prediction = list(fit.brm_exp3_causal_w,
                        fit.brm_exp3_causal_wh,
                        fit.brm_exp3_causal_whs) %>%
  map_dfr(~ fitted(., df.exp3.causal.means %>% 
                     left_join(df.exp3.model.causal,
                               by = c("clip", "ball")),
                   re_formula = NA) %>% 
            as_tibble()) %>% 
  mutate(ball = rep(c("A", "B"), n()/2),
         clip = rep(rep(1:32, each = 2), 3),
         model = rep(c("w", "wh", "whs"),
                     each = 64))

df.plot = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_both),
            by = c("clip")) %>% 
  left_join(model_prediction %>% 
              filter(model == model_index) %>% 
              select(model = Estimate, clip, ball),
            by = c("clip", "ball")) %>% 
  pivot_longer(cols = c(rating, model),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" &
                               outcome_both == 0, 1, colorindex),
         colorindex = ifelse(model == "rating" &
                               outcome_both == 1, 2, colorindex),
         colorindex = ifelse(model == "model", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "model"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>% 
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

p = actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
print(p)

# ggsave("../../figures/plots/exp3_causal_bars.pdf",
#        width = 12,
#        height = 6)

4.5.2.2 Selection

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
  

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(3, 7, 15, 16, 23)) %>%
  group_by(clip, ball) %>%
  summarize(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]) %>%
  left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
            by = c("clip", "ball")) %>%
  ungroup() %>%
  pivot_longer(cols = c(mean, w, wh, whs),
               names_to = "index",
               values_to = "value") %>% 
  mutate_at(.vars = vars(low, high),
            .funs = ~ ifelse(index == "mean", ., NA)) %>%
  mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
         clip = factor(clip, levels = c(7, 23, 3, 15, 16),
                       labels = c("causal chain",
                                  "double prevention",
                                  "joint causation",
                                  "overdetermination",
                                  "preemption")))

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         value = -10,
         index = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         index = NA,
         value = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = value,
                     group = index,
                     fill = index)) +
  geom_bar(stat = "identity", color = "black",
           position = position_dodge(0.9),
           width = 0.9) +
  geom_errorbar(mapping = aes(ymin = low, ymax = high),
                width = 0,
                size = 1,
                position = position_dodge(0.9)) +
  annotation_custom(grob = ball_a,
                    xmin = 0.5,
                    xmax = 1.5,
                    ymin = -25,
                    ymax = -7) + 
  annotation_custom(grob = ball_b,
                    xmin = 1.5,
                    xmax = 2.5,
                    ymin = -25,
                    ymax = -7) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) + 
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_grid(~clip) +
  scale_fill_grey(start = 1,
                  end = 0,
                  labels = c("mean rating  ", expression(CSM[W] ~ " "),
                             expression(CSM[WH] ~ " "),
                             expression(CSM[WHS] ~ " "))) +
  labs(y = "causal responsibility", fill = "") +
  coord_cartesian(clip = "off") + 
  theme_bw() +
  theme(legend.text.align = 0,
        text = element_text(size = 20),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 30),
        strip.text = element_text(size = 30),
        legend.spacing = unit(0.5, "cm"),
        legend.background = element_rect(fill = "transparent"),
        legend.margin = margin(t = 5, unit = "cm"),
        plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
print(p)

# ggsave("../../figures/plots/exp3_selection_bars.pdf",
#        width = 20,
#        height = 10,
#        device = cairo_pdf)

4.5.2.3 Scatter plot

func_scatterplot = function(model){
  if(model == "heuristic"){
    xlabel = "Heuristic"
  }else{
    xlabel = bquote(CSM[.(toupper(model))])
  }
  
  l.models = list(w = fit.brm_exp3_causal_w, 
                  wh = fit.brm_exp3_causal_wh,
                  whs = fit.brm_exp3_causal_whs,
                  heuristic = fit.brm_exp3_causal_heuristic)
  
  df.data = df.exp3.causal.means %>% 
    left_join(df.exp3.model.causal, by = c("clip", "ball")) %>% 
    left_join(df.regression.features, by = c("clip", "ball"))
  
  df.plot = l.models[[model]] %>%
    fitted(newdata = df.data,
           re_formula = NA) %>% 
    as_tibble() %>% 
    clean_names() %>%
    bind_cols(df.data)
  
  p = ggplot(data = df.plot,
             mapping = aes(x = estimate, 
                           y = rating_mean)) +
    geom_abline(intercept = 0,
                slope = 1,
                linetype = 2) + 
    geom_smooth(aes(y = estimate,
                    ymin = q2_5,
                    ymax = q97_5),
                stat = "identity",
                color = "black") +
    geom_linerange(size = 0.75, 
                   mapping = aes(ymin = rating_low,
                                 ymax = rating_high),
                   color = "gray50") +
    geom_point(size = 2) + 
    scale_color_manual(values = c("black", "#e41a1c"),
                       guide = F) +
    coord_fixed(xlim = c(0, 100),
                ylim = c(0, 100)) +
    labs(x = xlabel,
         y = "responsibility judgment") +
    annotate(geom = "text",
             label = str_c(
               "r = ", cor(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2) %>% 
                 as.character() %>% 
                 str_sub(start = 2),
               "\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2)),
             x = 5,
             y = 95,
             size = 6,
             hjust = 0) +
    scale_x_continuous(breaks = seq(0, 100, 25)) +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    theme(text = element_text(size = 20))
  return(p)
}

# for creating and saving an individual scatter plot 
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"

# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
#        width = 5,
#        height = 5)

list(func_scatterplot(model = "w"),
     func_scatterplot(model = "wh"),
     func_scatterplot(model = "whs"),
     func_scatterplot(model = "heuristic")) %>%
  wrap_plots(ncol = 2)

4.5.2.4 Individual participant variance for a selection of clips (clustered)

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(7, 23, 3, 15, 16)) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
    labels = c("causal chain", "double prevention", "joint causation",
      "overdetermination", "preemption")))

df.plot = df.plot %>% 
  left_join(df.cluster_whs %>% 
              group_by(cluster) %>% 
              mutate(group = factor(cluster,
                                    levels = 1:2,
                                    labels = str_c("n = ", group_size(.)))) %>% 
              ungroup() %>% 
              select(participant, cluster, group),
            by = "participant")


df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         rating = -10,
         group = NA,
         participant = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         group = NA,
         participant = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = rating,
                     group = participant,
                     shape = group)) +
  geom_line(mapping = aes(linetype = "individual",
                          color = group),
            size = 1,
            alpha = 0.3) +
  geom_point(mapping = aes(color = group),
             size = 1,
             alpha = 0.3) +
  stat_summary(mapping = aes(group = group,
                             color = group,
                             linetype = "mean"),
               fun = "mean",
               geom = "line",
               size = 1.5) +
  stat_summary(data = df.plot %>% filter(cluster == 1),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  stat_summary(data = df.plot %>% filter(cluster == 2),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  annotation_custom(grob = ball_a,
                    xmin = 0.5,
                    xmax = 1.5,
                    ymin = -30,
                    ymax = -10) + 
  annotation_custom(grob = ball_b,
                    xmin = 1.5,
                    xmax = 2.5,
                    ymin = -30,
                    ymax = -10) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) +
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_wrap(~ clip,
             ncol = 8) +
  coord_cartesian(xlim = c(0.9, 2.1),
                  ylim = c(-5, 105),
                  expand = F,
                  clip = "off") +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25)) +
  scale_color_brewer(palette = "Set1",
                     guide = "none") +
  labs(y = "causal responsibility rating",
       linetype = "legend",
       shape = "") +
  theme_bw() +
  guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
                                 keywidth = unit(1.2, "cm")),
         shape = guide_legend(override.aes = list(color = c("#377eb8",
                                                            "#e41a1c"),
                                                  shape = c(19, 19),
                                                  alpha = c(1, 1),
                                                  size = c(5, 5)))) +
  theme(panel.grid = element_blank(),
        text = element_text(size = 20),
        legend.position = c(0.505, 0.23),
        axis.title.x = element_blank(),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 24),
        legend.background = element_rect(fill = "transparent"),
        strip.text = element_text(size = 30),
        axis.text.x = element_blank(),
        legend.key = element_rect(fill = "transparent"),
        legend.box = "horizontal",
        legend.spacing.x = unit(0.1, "cm"),
        panel.spacing.x = unit(1, "cm"),
        plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))

print(p)

# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
#        plot = p,
#        width = 20,
#        height = 8.5,
#        device = cairo_pdf)

4.5.3 Tables

4.5.3.1 Relationship between predictors

df.exp3.model %>% 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(difference, whether, how, sufficient, robust) %>% 
  correlate() %>% 
  shave() %>% 
  fashion() %>% 
  print_table()
rowname difference whether how sufficient robust
difference
whether .50
how .79 .27
sufficient .21 .10 .36
robust .43 .93 .24 .20

4.5.3.2 Population level predictors in the mixed effects models

fit.brm_exp3_causal_w %>% 
  tidy(effects = "fixed") %>% 
  mutate(model = "CSM_w") %>% 
  bind_rows(fit.brm_exp3_causal_wh %>% 
              tidy(effects = "fixed") %>% 
              mutate(model = "CSM_wh")) %>% 
  bind_rows(fit.brm_exp3_causal_whs %>% 
              tidy(effects = "fixed") %>% 
              mutate(model = "CSM_whs")) %>% 
  mutate(term = tolower(term)) %>% 
  rename(`lower 95% HDI` = conf.low,
         `upper 95% HDI` = conf.high) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  select(model, everything(), -effect, -component) %>% 
  print_table()
model term estimate std.error lower 95% HDI upper 95% HDI
CSM_w (intercept) 31.68 1.61 28.41 34.87
CSM_w whether 39.30 2.16 35.11 43.54
CSM_wh (intercept) 10.29 1.56 7.30 13.41
CSM_wh whether 28.40 2.68 23.18 33.56
CSM_wh how 36.07 2.60 31.01 41.17
CSM_whs (intercept) 10.43 1.59 7.30 13.50
CSM_whs whether 27.93 2.60 22.80 33.09
CSM_whs how 26.47 2.91 20.83 32.42
CSM_whs sufficient 30.90 3.06 25.03 36.93

4.5.3.3 Individual participants regression results

df.exp3.causal.individual_fit %>% 
  select_if(~ (is.numeric(.) | is.factor(.))) %>% 
  mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>% 
  select(participant, everything(), best_model) %>% 
  print_table()
participant r_w r_wh r_whs rmse_w rmse_wh rmse_whs best_model
1 0.29 0.38 0.48 30.00 29.03 27.74 whs
2 0.20 0.77 0.77 39.32 27.94 27.55 whs
3 0.37 0.78 0.79 38.90 28.00 26.95 whs
4 0.41 0.55 0.63 43.62 40.60 38.31 whs
5 0.45 0.53 0.54 39.23 37.32 36.81 whs
6 0.18 0.54 0.54 40.84 36.32 36.07 whs
7 0.49 0.61 0.64 32.96 29.93 29.15 whs
8 0.39 0.63 0.68 38.37 33.26 31.64 whs
9 0.40 0.72 0.76 35.86 28.34 26.39 whs
10 0.28 0.52 0.53 29.13 26.28 25.85 whs
11 0.51 0.56 0.60 33.66 32.19 31.17 whs
12 0.56 0.57 0.70 32.78 32.08 28.84 whs
13 0.56 0.71 0.74 29.83 25.40 24.21 whs
14 0.43 0.47 0.50 33.72 32.82 32.23 whs
15 0.19 0.35 0.37 39.77 38.29 37.94 whs
16 0.62 0.67 0.76 40.49 37.98 34.39 whs
17 0.42 0.66 0.71 28.57 23.82 22.42 whs
18 0.31 0.35 0.39 37.78 37.09 36.57 whs
19 0.37 0.60 0.62 33.41 29.23 28.62 whs
20 0.57 0.61 0.68 36.47 35.00 32.71 whs
21 0.43 0.50 0.57 30.95 29.61 28.32 whs
22 0.51 0.67 0.71 33.72 29.36 28.00 whs
23 0.33 0.45 0.52 38.76 37.00 35.57 whs
24 0.46 0.63 0.69 34.01 29.96 28.22 whs
25 0.61 0.66 0.70 30.60 28.91 27.55 whs
26 0.57 0.71 0.76 40.43 35.03 32.80 whs
27 0.46 0.52 0.66 43.12 41.49 38.12 whs
28 0.30 0.45 0.46 29.33 27.62 27.38 whs
29 0.48 0.60 0.65 38.93 36.03 34.48 whs
30 0.18 0.70 0.70 33.06 24.95 24.71 whs
31 0.48 0.59 0.69 37.54 34.72 31.76 whs
32 0.37 0.61 0.70 40.62 35.94 32.92 whs
33 0.50 0.56 0.64 38.64 36.79 34.85 whs
34 0.24 0.51 0.52 44.50 40.86 40.19 whs
35 0.15 0.32 0.30 34.56 33.36 33.37 wh
36 0.19 0.84 0.82 45.38 29.58 29.70 wh
37 0.55 0.63 0.67 36.04 33.36 32.11 whs
38 0.52 0.61 0.67 36.03 33.39 31.53 whs
39 0.47 0.77 0.77 29.61 21.75 21.58 whs
40 0.46 0.78 0.79 32.41 23.60 22.69 whs
41 0.32 0.50 0.58 32.90 30.38 28.69 whs

4.5.3.4 CSMwhs predictions for selection of cases

df.exp3.model.causal %>%
  left_join(df.exp3.causal.regression,
            by = c("clip", "ball")) %>% 
  filter(clip %in% c(7, 23, 16, 3, 15)) %>%
  mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic),
            .funs = ~ round(., 2)) %>%
  select(clip, ball, difference, whether, how, sufficient, robust) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ ifelse(. < 0.5,
                             str_c("xmark (", ., ")"),
                             str_c("cmark (", ., ")"))) %>%
  arrange(clip) %>% 
  print_table()
clip ball difference whether how sufficient robust
7 A cmark (1) xmark (0.34) cmark (1) xmark (0) xmark (0.25)
7 B cmark (1) cmark (1) cmark (1) cmark (0.67) cmark (0.6)
23 A xmark (0.05) xmark (0) xmark (0) xmark (0) xmark (0)
23 B cmark (0.91) cmark (0.79) xmark (0) xmark (0) cmark (0.72)
16 A cmark (1) xmark (0.23) cmark (1) cmark (1) xmark (0.35)
16 B xmark (0) xmark (0) xmark (0) xmark (0) xmark (0)
3 A cmark (1) cmark (0.88) cmark (1) xmark (0.12) cmark (0.76)
3 B cmark (1) cmark (0.89) cmark (1) xmark (0.11) cmark (0.75)
15 A cmark (1) xmark (0.01) cmark (1) cmark (0.99) xmark (0.1)
15 B cmark (1) xmark (0.01) cmark (1) cmark (1) xmark (0.1)

4.5.3.5 CSMwhs predictions for all cases

df.exp3.causal.regression %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball")) %>% 
  left_join(df.exp3.clipinfo %>% 
              select(-clipindex),
            by = c("clip")) %>% 
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ . * 100) %>% 
  mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
  select(clip, ball, contains("outcome"), difference, whether, how, sufficient,
         robust, w, wh, whs, heuristic, rating = rating_mean) %>% 
  # print_table(format = "latex", digits = 0)
  print_table(digits = 0)
clip ball outcome_both outcome_a outcome_b outcome_none difference whether how sufficient robust w wh whs heuristic rating
1 A 0 0 0 0 100 40 100 23 36 47 58 55 57 42
1 B 0 0 0 0 100 15 100 16 9 38 51 46 54 37
2 A 0 0 0 0 57 12 0 0 10 36 14 14 25 21
2 B 0 0 0 0 18 0 0 0 0 32 10 10 24 19
3 A 1 0 0 0 100 88 100 12 76 66 71 65 72 76
3 B 1 0 0 0 100 89 100 11 75 67 72 65 72 75
4 A 1 0 0 0 100 78 100 4 78 62 69 60 58 63
4 B 1 0 0 0 100 95 100 15 57 69 73 68 54 78
5 A 0 0 1 0 100 90 100 0 47 67 72 62 47 71
5 B 0 0 1 0 100 0 100 0 0 32 46 37 68 22
6 A 0 0 1 0 100 59 100 16 35 55 63 59 53 73
6 B 0 0 1 0 100 18 100 6 14 39 52 44 53 22
7 A 1 0 1 0 100 34 100 0 25 45 56 46 70 59
7 B 1 0 1 0 100 100 100 67 60 71 75 86 64 79
8 A 1 0 1 0 0 0 0 0 0 32 10 10 25 7
8 B 1 0 1 0 100 100 100 100 100 71 75 96 84 92
9 A 0 1 0 0 0 0 0 0 0 32 10 10 14 8
9 B 0 1 0 0 100 100 100 0 100 71 75 65 78 90
10 A 0 1 0 0 77 18 0 0 22 39 15 16 15 23
10 B 0 1 0 0 98 79 0 0 63 63 33 33 21 55
11 A 1 1 0 0 100 70 100 77 68 59 66 80 71 93
11 B 1 1 0 0 0 0 0 0 0 32 10 10 16 4
12 A 1 1 0 0 100 82 100 74 83 64 70 83 57 77
12 B 1 1 0 0 100 0 100 12 24 32 46 41 53 37
13 A 0 1 1 0 67 34 0 0 35 45 20 20 14 8
13 B 0 1 1 0 70 35 0 0 35 46 20 20 21 64
14 A 0 1 1 0 97 91 0 0 59 67 36 36 21 22
14 B 0 1 1 0 91 77 0 0 51 62 32 32 20 18
15 A 1 1 1 0 100 1 100 99 10 32 47 68 71 76
15 B 1 1 1 0 100 1 100 100 10 32 47 68 71 76
16 A 1 1 1 0 100 23 100 100 35 41 53 74 80 92
16 B 1 1 1 0 0 0 0 0 0 32 10 10 14 4
17 A 0 0 0 1 100 19 100 37 18 39 52 54 66 69
17 B 0 0 0 1 100 0 100 36 17 32 46 48 65 46
18 A 0 0 0 1 100 11 100 40 17 36 49 52 55 63
18 B 0 0 0 1 100 7 100 37 9 34 48 50 56 66
19 A 1 0 0 1 100 74 100 7 65 61 67 60 55 53
19 B 1 0 0 1 100 72 100 7 65 60 67 59 55 49
20 A 1 0 0 1 100 92 100 8 72 68 72 65 57 41
20 B 1 0 0 1 100 88 100 4 53 66 71 63 56 71
21 A 0 0 1 1 100 47 100 40 45 50 60 63 58 80
21 B 0 0 1 1 100 9 100 21 10 35 49 46 59 18
22 A 0 0 1 1 100 100 100 89 83 71 75 92 48 60
22 B 0 0 1 1 100 8 100 0 15 35 49 39 53 42
23 A 1 0 1 1 5 0 0 0 0 32 10 10 15 3
23 B 1 0 1 1 91 79 0 0 72 63 33 33 22 39
24 A 1 0 1 1 100 66 100 4 63 58 65 57 57 44
24 B 1 0 1 1 100 94 100 22 79 69 73 70 54 73
25 A 0 1 0 1 100 25 100 21 26 42 53 50 69 43
25 B 0 1 0 1 100 74 100 54 65 61 67 74 56 73
26 A 0 1 0 1 100 6 100 3 9 34 48 40 60 39
26 B 0 1 0 1 100 87 100 35 54 66 71 72 46 69
27 A 1 1 0 1 100 97 100 52 97 70 74 80 67 80
27 B 1 1 0 1 0 0 0 0 0 32 10 10 17 6
28 A 1 1 0 1 100 90 100 22 80 67 72 69 79 89
28 B 1 1 0 1 0 0 0 0 0 32 10 10 12 5
29 A 0 1 1 1 100 58 100 24 44 55 63 61 66 47
29 B 0 1 1 1 100 63 100 24 38 57 64 62 54 67
30 A 0 1 1 1 100 57 100 29 49 54 63 62 63 58
30 B 0 1 1 1 100 46 100 24 39 50 60 57 63 56
31 A 1 1 1 1 100 2 100 4 3 32 47 39 63 44
31 B 1 1 1 1 100 4 100 4 4 33 47 39 53 46
32 A 1 1 1 1 0 0 0 0 0 32 10 10 16 5
32 B 1 1 1 1 100 75 100 66 73 61 68 78 65 71

4.5.3.6 Heuristic model

fit.brm_exp3_causal_heuristic %>% 
  tidy(effects = "fixed") %>% 
  mutate(term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = conf.low,
         `upper 95% HDI` = conf.high) %>% 
  print_table()
Warning in tidy.brmsfit(., effects = "fixed"): some parameter names contain
underscores: term naming may be unreliable!
effect component term estimate std.error lower 95% HDI upper 95% HDI
fixed cond (intercept) 49.79 1.46 46.91 52.74
fixed cond moving 0.22 0.22 0.01 0.81
fixed cond speed 2.06 0.85 0.44 3.67
fixed cond contact_e 0.38 0.35 0.01 1.29
fixed cond e_speed_diff 0.12 0.12 0.00 0.46
fixed cond e_direction_diff 1.04 0.74 0.04 2.71
fixed cond total_speed_diff 2.21 0.95 0.41 4.17
fixed cond total_direction_diff 3.99 0.91 2.14 5.68
fixed cond transfer 15.59 0.79 14.06 17.14
fixed cond e_moving 0.18 0.17 0.01 0.62
fixed cond exclusive 4.39 0.73 2.97 5.81

5 Appendix

5.1 Experiment 3

5.1.1 Does the question framing matter? (Responsibility vs. causation)

We have run a pilot eye-tracking experiment with 15 participants who saw the same clips as the ones in Experiment 3. Instead of asking participants to judge “To what extent were A and B responsible for E (not) going through the gate?”, we asked them to judge to what extent they agreed with the following statement: “Ball A/B caused ball E to go through the gate.” or “Ball A/B prevented ball E from going through the gate.” depending on the outcome.

Participants’ agreement judgments with the causal statements in this pilot eye-tracking experiment were strikingly similar to participants’ responsibility judgments in the online experiment.

df.exp3.eye_tracking_pilot = 
  read_csv("../../data/experiment3_eye_tracking_pilot.csv") %>% 
  rename(clip = trial) %>%
  mutate(ball = str_to_upper(ball)) %>% 
  group_by(clip, ball, outcome) %>% 
  summarize(causality_mean = mean(causation),
            causality_low = smean.cl.boot(causation)[2],
            causality_high = smean.cl.boot(causation)[3]) %>% 
  ungroup()

df.plot = df.exp3.eye_tracking_pilot %>% 
  # select(clip, ball, outcome,) %>% 
  left_join(df.exp3.causal.means %>% 
              select(clip,
                     ball,
                     responsibility_mean = rating_mean,
                     responsibility_low = rating_low,
                     responsibility_high = rating_high),
            by = c("clip", "ball"))

# highlight the double prevention clip (#23)
df.plot = df.plot %>% 
  mutate(color = ifelse(clip == 23, "1", "0"))

ggplot(data = df.plot,
       mapping = aes(x = causality_mean,
                     y = responsibility_mean)) + 
  geom_abline(intercept = 0, 
              slope = 1,
              linetype = 2) + 
  geom_linerange(mapping = aes(ymin = responsibility_low,
                               ymax = responsibility_high),
                 alpha = 0.1) +
  geom_linerange(mapping = aes(xmin = causality_low,
                               xmax = causality_high),
                 alpha = 0.1) +
  geom_point(size = 2,
             mapping = aes(color = color),
             show.legend = F) + 
  geom_smooth(method = "lm",
              color = "black") +
  annotate(geom = "text",
           x = c(0, 0),
           y = c(100, 90),
           label = c(str_c("r = ",
                           cor(df.plot$causality_mean,
                               df.plot$responsibility_mean) %>% 
                             round(2)),
                     str_c("RMSE = ", rmse(df.plot$causality_mean,
                                           df.plot$responsibility_mean) %>% 
                             round(2))),
           hjust = 0, 
           size = 8) + 
  scale_x_continuous(breaks = seq(0, 100, 25),
                     limits = c(0, 100)) +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     limits = c(0, 100)) +
  scale_color_manual(values = c("black", "red")) +
  labs(x = "causal judgment",
       y = "responsibility judgment") +
  theme(text = element_text(size = 24))

# ggsave("../../figures/plots/aux_question_framing_comparison.pdf",
#        width = 6,
#        height = 6)

5.1.2 CSM_WHS fits for each individual participant

df.plot = fit.brm_exp3_causal_whs %>% 
  fitted() %>% 
  as_tibble() %>% 
  select(prediction = Estimate) %>% 
  bind_cols(df.exp3.causal.long) %>% 
  relocate(prediction, .after = last_col())

df.text = df.plot %>% 
  group_by(participant) %>% 
  summarize(r = round(cor(prediction, rating), 2)) %>% 
  ungroup() %>% 
  mutate(r = str_c("r = ", r),
         prediction = 1,
         rating = 110)

ggplot(data = df.plot,
       mapping = aes(x = prediction,
                     y = rating)) + 
  geom_abline(intercept = 0, 
              slope = 1,
              color = "blue",
              alpha = 0.5) +
  geom_point(alpha = 0.3,
             size = 1) +
  geom_text(data = df.text,
            mapping = aes(label = r),
            hjust = 0,
            color = "red",
            size = 4) +
  facet_wrap(facets = vars(participant),
             ncol = 7) +
  labs(x = "model prediction",
       y = "causal responsibility rating") + 
  coord_cartesian(xlim = c(0, 100),
                  ylim = c(0, 100),
                  expand = F,
                  clip = "off") +
  scale_x_continuous(breaks = seq(0, 100, 25),
                     limits = c(0, 100)) +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     expand = expansion(mult = c(0, 0))) +
  theme_classic() + 
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        text = element_text(size = 16),
        panel.spacing.x = unit(0.75, "cm"),
        panel.spacing.y = unit(1, "cm"),
        plot.margin = margin(t = 0.7,
                             r = 0.8,
                             b = 0.2,
                             l = 0.2,
                             unit = "cm"))

# ggsave("../../figures/plots/exp3_individual_scatter_plots.pdf",
#        width = 12,
#        height = 8)

5.1.3 Individual participant fit model comparison

df.plot = df.exp3.causal.individual_fit %>% 
  select(participant, contains("r_"), contains("rmse_")) %>% 
  pivot_longer(cols = -participant,
               names_to = c("measure", "model"),
               names_sep = "_",
               values_to = "fit")

ggplot(data = df.plot %>% 
         filter(measure == "r"),
       mapping = aes(x = fit,
                     fill = model)) + 
  geom_density(alpha = 0.5) + 
  labs(x = "correlation value") + 
  scale_x_continuous(breaks = seq(0, 1, 0.2),
                     limits = c(0, 1)) + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) + 
  scale_fill_brewer(palette = "Set1") + 
  # ggplot2::theme_classic() + 
  theme(legend.position = c(0.3, 0.95),
        legend.direction = "horizontal",
        text = element_text(size = 20))

# ggsave("../../figures/plots/exp3_individual_densities.pdf",
#        width = 8,
#        height = 6)

### 3D scatter plot of participant clusters

plot_ly(x = df.cluster_whs$whether,
        y = df.cluster_whs$how,
        z = df.cluster_whs$sufficient,
        type = "scatter3d",
        mode = "markers",
        color = as.factor(df.cluster_whs$cluster),
        colors = c("#e41a1c", "#377eb8")) %>% 
  layout(showlegend = FALSE,
         title = "Participant clusters",
         scene = list(
           xaxis = list(title = "whether"),
           yaxis = list(title = "how"),
           zaxis = list(title = "sufficient")))

5.1.4 Clusters of participants

set.seed(1)

df.plot = df.exp3.causal.long %>% 
  left_join(df.cluster_whs %>% 
              select(participant, cluster) %>% 
              mutate(participant = as.numeric(participant),
                     cluster = as.factor(cluster)) %>% 
              group_by(cluster) %>% 
              mutate(label = factor(cluster,
                                    levels = 1:2,
                                    labels = str_c("n = ", group_size(.)))) %>% 
              ungroup(),
            by = "participant")

ggplot(data = df.plot,
       mapping = aes(x = ball,
                     y = rating,
                     group = label,
                     fill = label)) + 
  geom_point(shape = 21,
             alpha = 0.2,
             position = position_jitterdodge(jitter.width = 0.1,
                                             jitter.height = 0,
                                             dodge.width = 0.75)) + 
  stat_summary(fun.data = "mean_cl_boot",
               geom = "linerange", 
               mapping = aes(color = label),
               size = 1,
               position = position_dodge(width = 0.75)) + 
  stat_summary(fun = "mean",
               geom = "point", 
               size = 4,
               shape = 21,
               position = position_dodge(width = 0.75)) + 
  facet_wrap(facets = vars(clip), ncol = 8) + 
  labs(y = "causal responsibility rating",
       fill = "cluster",
       color = "cluster") + 
  scale_fill_brewer(palette = "Set1") + 
  scale_color_brewer(palette = "Set1") + 
  theme(text = element_text(size = 24),
        legend.position = "bottom")

# ggsave("../../figures/plots/exp3_clusters_points.pdf",
#        width = 12,
#        height = 8)

5.1.5 Ternary plots for individual regression results

library("ggtern")

df.plot = df.exp3.causal.individual_fit %>% 
  select(participant, fit_whs) %>% 
  mutate(estimates = map(fit_whs, ~ fixef(.) %>% 
                           as_tibble(rownames = "term") %>% 
                           clean_names())) %>%
  select(participant, estimates) %>% 
  unnest(estimates) %>% 
  filter(!str_detect(term, "Intercept")) %>% 
  select(participant, term, estimate) %>% 
  pivot_wider(names_from = term,
              values_from = estimate) %>% 
  mutate_at(vars(whether, how, sufficient),
            .funs = list(norm = ~ . / (how + whether + sufficient))) %>% 
  mutate(color = 0,
         color = ifelse(how_norm == max(how_norm), 1, color),
         color = ifelse(whether_norm == max(whether_norm), 2, color),
         color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
         color = factor(color))

ggplot(data = df.plot,
       mapping = aes(x = how,
                     y = sufficient,
                     z = whether,
                     color = color)) + 
  geom_point(alpha = 0.7,
             size = 2,
             show.legend = F) + 
  scale_color_manual(values = c("black", "red", "green", "blue")) + 
  coord_tern() +
  theme_showarrows() + 
  theme(text = element_text(size = 20),
        tern.axis.title.T = element_blank(),
        tern.axis.title.L = element_blank(),
        tern.axis.title.R = element_blank())

# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
#        width = 5,
#        height = 5)

6 Session info

sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_0.5.0     stringr_1.4.0     dplyr_1.0.0       purrr_0.3.4      
 [5] readr_1.3.1       tidyr_1.1.0       tibble_3.0.1      tidyverse_1.3.0  
 [9] patchwork_1.0.0   clValid_0.6-7     cluster_2.1.0     egg_0.4.5        
[13] gridExtra_2.3     plotly_4.9.2.1    png_0.1-7         tidybayes_2.1.1  
[17] Hmisc_4.4-0       ggplot2_3.3.1     Formula_1.2-3     survival_3.1-8   
[21] lattice_0.20-38   janitor_2.0.1     jsonlite_1.6.1    xtable_1.8-4     
[25] corrr_0.4.2       brms_2.13.0       Rcpp_1.0.4.6      broom.mixed_0.2.6
[29] lme4_1.1-23       Matrix_1.2-18     kableExtra_1.1.0  knitr_1.28       

loaded via a namespace (and not attached):
  [1] readxl_1.3.1         backports_1.1.7      plyr_1.8.6          
  [4] igraph_1.2.5         lazyeval_0.2.2       svUnit_1.0.3        
  [7] TMB_1.7.16           splines_3.6.3        crosstalk_1.1.0.1   
 [10] rstantools_2.1.0     inline_0.3.15        digest_0.6.25       
 [13] htmltools_0.4.0      rsconnect_0.8.16     fansi_0.4.1         
 [16] magrittr_1.5         checkmate_2.0.0      modelr_0.1.8        
 [19] RcppParallel_5.0.1   matrixStats_0.56.0   xts_0.12-0          
 [22] prettyunits_1.1.1    jpeg_0.1-8.1         colorspace_1.4-1    
 [25] blob_1.2.1           rvest_0.3.5          ggdist_2.1.1        
 [28] haven_2.3.1          xfun_0.14            callr_3.4.3         
 [31] crayon_1.3.4         zoo_1.8-8            glue_1.4.1          
 [34] gtable_0.3.0         emmeans_1.4.7        webshot_0.5.2       
 [37] pkgbuild_1.0.8       rstan_2.19.3         abind_1.4-5         
 [40] scales_1.1.1         mvtnorm_1.1-1        DBI_1.1.0           
 [43] miniUI_0.1.1.1       viridisLite_0.3.0    htmlTable_1.13.3    
 [46] HDInterval_0.2.2     foreign_0.8-75       stats4_3.6.3        
 [49] StanHeaders_2.21.0-5 DT_0.13              htmlwidgets_1.5.1   
 [52] httr_1.4.1           threejs_0.3.3        arrayhelpers_1.1-0  
 [55] RColorBrewer_1.1-2   acepack_1.4.1        ellipsis_0.3.1      
 [58] farver_2.0.3         pkgconfig_2.0.3      loo_2.2.0           
 [61] dbplyr_1.4.4         nnet_7.3-12          utf8_1.1.4          
 [64] labeling_0.3         tidyselect_1.1.0     rlang_0.4.6         
 [67] reshape2_1.4.4       later_1.1.0.1        cellranger_1.1.0    
 [70] munsell_0.5.0        tools_3.6.3          cli_2.0.2           
 [73] generics_0.0.2       broom_0.7.0.9000     ggridges_0.5.2      
 [76] evaluate_0.14        fastmap_1.0.1        yaml_2.2.1          
 [79] fs_1.4.1             processx_3.4.2       nlme_3.1-144        
 [82] mime_0.9             xml2_1.3.2           compiler_3.6.3      
 [85] bayesplot_1.7.2      shinythemes_1.1.2    rstudioapi_0.11     
 [88] reprex_0.3.0         statmod_1.4.34       stringi_1.4.6       
 [91] highr_0.8            ps_1.3.3             Brobdingnag_1.2-6   
 [94] nloptr_1.2.2.1       markdown_1.1         shinyjs_1.1         
 [97] vctrs_0.3.1          pillar_1.4.4         lifecycle_0.2.0     
[100] bridgesampling_1.0-0 estimability_1.3     data.table_1.12.8   
[103] httpuv_1.5.4         R6_2.4.1             latticeExtra_0.6-29 
[106] bookdown_0.19        promises_1.1.1       boot_1.3-24         
[109] colourpicker_1.0     MASS_7.3-51.5        gtools_3.8.2        
[112] assertthat_0.2.1     withr_2.2.0          shinystan_2.5.0     
[115] mgcv_1.8-31          parallel_3.6.3       hms_0.5.3           
[118] rpart_4.1-15         class_7.3-15         coda_0.19-3         
[121] minqa_1.2.4          rmarkdown_2.2        snakecase_0.11.0    
[124] shiny_1.4.0.2        lubridate_1.7.9      base64enc_0.1-3     
[127] dygraphs_1.1.1.6